Nrf2 KO (-/-) PEITC 16S Microbiome Data Visualization

# NOTE (DS, 02/09/2018): mean plots were incorrect; using new function
require(data.table)
require(phyloseq)
require(ggplot2)
require(plotly)
require(DT)
require(shiny)
# source("source/functions_v2.R")
source("source/functions_may2019.R")

Load data

# Counts
load("data/ps.RData")
# Taxonomy
load("data/taxa.plus.RData")
taxa <- data.table(seq16s = rownames(taxa.plus),
                   taxa.plus)
# Samples
# ps@sam_data
load("data/samples.RData")
samples$Sample <- substr(x = samples$Name,
                         start = 1,
                         stop = 5)
samples$Sample[samples$Sample %in% c("4A_S1",
                                     "4B_S2",
                                     "4C_S3")] <- 
  substr(x = samples$Sample[samples$Sample %in% c("4A_S1",
                                                  "4B_S2",
                                                  "4C_S3")],
         start = 1,
         stop = 2)
DT::datatable(samples)

Taxonomy Ranks:

King Phillip Can nOt Find Green Socks

• Kingdom
  
• Phylum
  
• Class
  
• Order
  
• Family
  
• Genus
  
• Species

Richness

# Part I: keep controls and PEITCs at weeks 5 and 9 only
w59 <- prune_samples(samples = sample_names(ps)[!(sample_names(ps) %in%
                                                    c("4A",
                                                      "4B",
                                                      "4C",
                                                      "Undetermined"))],
                    x = ps)
# w59 <- prune_samples(samples = sample_names(ps)[sample_names(ps) != "Undetermined"], 
#                     x = ps)
# p1 <- plot_richness(w59,
#               x = "Diet_Week", 
#               measures = "Shannon",
#               color = "Week") +
#   geom_line(aes(group = ID),
#             color = "black") +
#   geom_point(shape = 21,
#              size = 3,
#              color = "black")
p1 <- plot_richness(w59,
                    x = "Diet_Week", 
                    measures = "Shannon") +
    geom_line(aes(group = ID),
            color = "black") +
  geom_point(aes(fill = ID),
             shape = 21,
             size = 3,
             color = "black") +
  scale_x_discrete("") +
  theme(axis.text.x = element_text(angle = 30,
                                   hjust = 1,
                                   vjust = 1),
        legend.position = "none")
ggplotly(p1)

tiff(filename = "tmp/ko_shannon.tiff",
     height = 4,
     width = 5,
     units = "in",
     res = 600,
     compression = "lzw+p")
print(p1)
graphics.off()

OTU table (first 10 rows)

Total counts per sample (i.e. sequencing depth)

Counts at Phylum level

Relative abundance (%) at Phylum level

Remove phyla with relative abundance of >= 1% in less than 10% of samples.

t1 <- data.table(Phylum = ra_p$Phylum,
                 `Number of Samples` = rowSums(ra_p[, 2:ncol(ra_p)] >= 0.01))
t1$`Percent Samples` <-  t1$`Number of Samples`/30
setorder(t1, -`Number of Samples`)
datatable(t1,
          rownames = FALSE,
          caption = "Taxonomic  count table",
          class = "cell-border stripe",
          options = list(search = FALSE,
                         pageLength = nrow(t1))) %>%
  formatPercentage(columns = 3,
                   digits = 1)

[1] "Bacteroidetes, Firmicutes, Proteobacteria, Actinobacteria, Saccharibacteria, Deferribacteres, Cyanobacteria"

Relative Abundance in Samples at Different Taxonomic Ranks

1. Class

p0 <- ggplot(mu,
             aes(x = Week,
                 y = x,
                 group = Treatment)) +
  facet_wrap(~ Class,
             scale = "free_y") +
  geom_line() +
  geom_point(aes(shape = Treatment,
                 color = Treatment),
             size = 3,
             alpha = 0.5) +
  scale_x_discrete("") +
  scale_y_continuous("Relative Abundance (%)") +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 45,
                                   hjust = 1))
tiff(filename = "tmp/ko_class_over_time.tiff",
     height = 5,
     width = 7,
     units = "in",
     res = 600,
     compression = "lzw+p")
print(p0)
graphics.off()
print(p0)

p1 <- ggplot(mu,
             aes(x = x,
                 y = Class,
                 color = Treatment,
                 shape = Week)) +
  geom_point(size = 3,
             alpha = 0.5) +
  geom_vline(xintercept = 1,
             linetype = "dashed") +
  scale_x_continuous("Relative Abundance (%)") +
  theme(legend.position = "top")
tiff(filename = "tmp/ko_class_ra.tiff",
     height = 4,
     width = 7,
     units = "in",
     res = 600,
     compression = "lzw+p")
print(p1)
graphics.off()
ggplotly(p1)